home *** CD-ROM | disk | FTP | other *** search
/ The World of Computer Software / The World of Computer Software.iso / srcuc.zip / BIGNUM.C < prev    next >
C/C++ Source or Header  |  1991-10-29  |  50KB  |  1,674 lines
  1. /* -*-C-*-
  2.  
  3. $Header: /scheme/users/jinx/microcode/RCS/bignum.c,v 9.36 1991/10/29 22:55:11 jinx Exp $
  4.  
  5. Copyright (c) 1989-91 Massachusetts Institute of Technology
  6.  
  7. This material was developed by the Scheme project at the Massachusetts
  8. Institute of Technology, Department of Electrical Engineering and
  9. Computer Science.  Permission to copy this software, to redistribute
  10. it, and to use it for any purpose is granted, subject to the following
  11. restrictions and understandings.
  12.  
  13. 1. Any copy made of this software must include this copyright notice
  14. in full.
  15.  
  16. 2. Users of this software agree to make their best efforts (a) to
  17. return to the MIT Scheme project any improvements or extensions that
  18. they make, so that these may be included in future releases; and (b)
  19. to inform MIT of noteworthy uses of this software.
  20.  
  21. 3. All materials developed as a consequence of the use of this
  22. software shall duly acknowledge such use, in accordance with the usual
  23. standards of acknowledging credit in academic research.
  24.  
  25. 4. MIT has made no warrantee or representation that the operation of
  26. this software will be error-free, and MIT is under no obligation to
  27. provide any services, by way of maintenance, update, or otherwise.
  28.  
  29. 5. In conjunction with products arising from the use of this material,
  30. there shall be no use of the name of the Massachusetts Institute of
  31. Technology nor of any adaptation thereof in any advertising,
  32. promotional, or sales literature without prior written consent from
  33. MIT in each case. */
  34.  
  35. /* Implementation of Bignums (unlimited precision integers) */
  36.  
  37. #ifdef MIT_SCHEME
  38. #include "scheme.h"
  39. #undef CHAR_BIT /* redefined in "limits.h" */
  40. #else
  41. #include "bignum.h"
  42. #endif
  43. #include "bignumint.h"
  44. #include "limits.h"
  45.  
  46. #ifndef MIT_SCHEME
  47.  
  48. static bignum_type
  49. DEFUN (bignum_malloc, (length), bignum_length_type length)
  50. {
  51.   extern char * malloc ();
  52.   char * result = (malloc ((length + 1) * (sizeof (bignum_digit_type))));
  53.   BIGNUM_ASSERT (result != ((char *) 0));
  54.   return ((bignum_type) result);
  55. }
  56.  
  57. static bignum_type
  58. DEFUN (bignum_realloc, (bignum, length),
  59.        bignum_type bignum AND bignum_length_type length)
  60. {
  61.   extern char * realloc ();
  62.   char * result =
  63.     (realloc (((char *) bignum),
  64.           ((length + 1) * (sizeof (bignum_digit_type)))));
  65.   BIGNUM_ASSERT (result != ((char *) 0));
  66.   return ((bignum_type) result);
  67. }
  68.  
  69. #endif /* not MIT_SCHEME */
  70.  
  71. /* Forward references */
  72. static int EXFUN (bignum_equal_p_unsigned,
  73.           (bignum_type, bignum_type));
  74. static enum bignum_comparison EXFUN (bignum_compare_unsigned,
  75.                      (bignum_type, bignum_type));
  76. static bignum_type EXFUN (bignum_add_unsigned,
  77.               (bignum_type, bignum_type, int));
  78. static bignum_type EXFUN (bignum_subtract_unsigned,
  79.               (bignum_type, bignum_type));
  80. static bignum_type EXFUN (bignum_multiply_unsigned,
  81.               (bignum_type, bignum_type, int));
  82. static bignum_type EXFUN (bignum_multiply_unsigned_small_factor,
  83.               (bignum_type, bignum_digit_type, int));
  84. static void EXFUN (bignum_destructive_scale_up,
  85.            (bignum_type, bignum_digit_type));
  86. static void EXFUN (bignum_destructive_add,
  87.            (bignum_type, bignum_digit_type));
  88. static void EXFUN (bignum_divide_unsigned_large_denominator,
  89.            (bignum_type, bignum_type, bignum_type *, bignum_type *,
  90.             int, int));
  91. static void EXFUN (bignum_destructive_normalization,
  92.            (bignum_type, bignum_type, int));
  93. static void EXFUN (bignum_destructive_unnormalization,
  94.            (bignum_type, int));
  95. static void EXFUN (bignum_divide_unsigned_normalized,
  96.            (bignum_type, bignum_type, bignum_type));
  97. static bignum_digit_type EXFUN (bignum_divide_subtract,
  98.                 (bignum_digit_type *, bignum_digit_type *,
  99.                  bignum_digit_type, bignum_digit_type *));
  100. static void EXFUN (bignum_divide_unsigned_medium_denominator,
  101.            (bignum_type, bignum_digit_type, bignum_type *,
  102.             bignum_type *, int, int));
  103. static bignum_digit_type EXFUN (bignum_digit_divide,
  104.                 (bignum_digit_type, bignum_digit_type,
  105.                  bignum_digit_type, bignum_digit_type *));
  106. static bignum_digit_type EXFUN (bignum_digit_divide_subtract,
  107.                 (bignum_digit_type, bignum_digit_type,
  108.                  bignum_digit_type, bignum_digit_type *));
  109. static void EXFUN (bignum_divide_unsigned_small_denominator,
  110.            (bignum_type, bignum_digit_type, bignum_type *,
  111.             bignum_type *, int, int));            
  112. static bignum_digit_type EXFUN (bignum_destructive_scale_down,
  113.                 (bignum_type, bignum_digit_type));
  114. static bignum_type EXFUN (bignum_remainder_unsigned_small_denominator,
  115.               (bignum_type, bignum_digit_type, int));
  116. static bignum_type EXFUN (bignum_digit_to_bignum,
  117.               (bignum_digit_type, int));
  118. static bignum_type EXFUN (bignum_allocate,
  119.               (bignum_length_type, int));
  120. static bignum_type EXFUN (bignum_allocate_zeroed,
  121.               (bignum_length_type, int));
  122. static bignum_type EXFUN (bignum_shorten_length,
  123.               (bignum_type, bignum_length_type));
  124. static bignum_type EXFUN (bignum_trim,
  125.               (bignum_type));
  126. static bignum_type EXFUN (bignum_copy,
  127.               (bignum_type));
  128. static bignum_type EXFUN (bignum_new_sign,
  129.               (bignum_type, int));
  130. static bignum_type EXFUN (bignum_maybe_new_sign,
  131.               (bignum_type, int));
  132. static void EXFUN (bignum_destructive_copy,
  133.            (bignum_type, bignum_type));
  134. static void EXFUN (bignum_destructive_zero,
  135.            (bignum_type));
  136.  
  137. /* Exports */
  138.  
  139. bignum_type
  140. DEFUN_VOID (bignum_make_zero)
  141. {
  142.   fast bignum_type result = (BIGNUM_ALLOCATE (0));
  143.   BIGNUM_SET_HEADER (result, 0, 0);
  144.   return (result);
  145. }
  146.  
  147. bignum_type
  148. DEFUN (bignum_make_one, (negative_p), int negative_p)
  149. {
  150.   fast bignum_type result = (BIGNUM_ALLOCATE (1));
  151.   BIGNUM_SET_HEADER (result, 1, negative_p);
  152.   (BIGNUM_REF (result, 0)) = 1;
  153.   return (result);
  154. }
  155.  
  156. int
  157. DEFUN (bignum_equal_p, (x, y),
  158.        fast bignum_type x AND fast bignum_type y)
  159. {
  160.   return
  161.     ((BIGNUM_ZERO_P (x))
  162.      ? (BIGNUM_ZERO_P (y))
  163.      : ((! (BIGNUM_ZERO_P (y)))
  164.     && ((BIGNUM_NEGATIVE_P (x))
  165.         ? (BIGNUM_NEGATIVE_P (y))
  166.         : (! (BIGNUM_NEGATIVE_P (y))))
  167.     && (bignum_equal_p_unsigned (x, y))));
  168. }
  169.  
  170. enum bignum_comparison
  171. DEFUN (bignum_test, (bignum), fast bignum_type bignum)
  172. {
  173.   return
  174.     ((BIGNUM_ZERO_P (bignum))
  175.      ? bignum_comparison_equal
  176.      : (BIGNUM_NEGATIVE_P (bignum))
  177.      ? bignum_comparison_less
  178.      : bignum_comparison_greater);
  179. }
  180.  
  181. enum bignum_comparison
  182. DEFUN (bignum_compare, (x, y),
  183.        fast bignum_type x AND fast bignum_type y)
  184. {
  185.   return
  186.     ((BIGNUM_ZERO_P (x))
  187.      ? ((BIGNUM_ZERO_P (y))
  188.     ? bignum_comparison_equal
  189.     : (BIGNUM_NEGATIVE_P (y))
  190.     ? bignum_comparison_greater
  191.     : bignum_comparison_less)
  192.      : (BIGNUM_ZERO_P (y))
  193.      ? ((BIGNUM_NEGATIVE_P (x))
  194.     ? bignum_comparison_less
  195.     : bignum_comparison_greater)
  196.      : (BIGNUM_NEGATIVE_P (x))
  197.      ? ((BIGNUM_NEGATIVE_P (y))
  198.     ? (bignum_compare_unsigned (y, x))
  199.     : (bignum_comparison_less))
  200.      : ((BIGNUM_NEGATIVE_P (y))
  201.     ? (bignum_comparison_greater)
  202.     : (bignum_compare_unsigned (x, y))));
  203. }
  204.  
  205. bignum_type
  206. DEFUN (bignum_add, (x, y),
  207.        fast bignum_type x AND fast bignum_type y)
  208. {
  209.   return
  210.     ((BIGNUM_ZERO_P (x))
  211.      ? (BIGNUM_MAYBE_COPY (y))
  212.      : (BIGNUM_ZERO_P (y))
  213.      ? (BIGNUM_MAYBE_COPY (x))
  214.      : ((BIGNUM_NEGATIVE_P (x))
  215.     ? ((BIGNUM_NEGATIVE_P (y))
  216.        ? (bignum_add_unsigned (x, y, 1))
  217.        : (bignum_subtract_unsigned (y, x)))
  218.     : ((BIGNUM_NEGATIVE_P (y))
  219.        ? (bignum_subtract_unsigned (x, y))
  220.        : (bignum_add_unsigned (x, y, 0)))));
  221. }
  222.  
  223. bignum_type
  224. DEFUN (bignum_subtract, (x, y),
  225.        fast bignum_type x AND fast bignum_type y)
  226. {
  227.   return
  228.     ((BIGNUM_ZERO_P (x))
  229.      ? ((BIGNUM_ZERO_P (y))
  230.     ? (BIGNUM_MAYBE_COPY (y))
  231.     : (bignum_new_sign (y, (! (BIGNUM_NEGATIVE_P (y))))))
  232.      : ((BIGNUM_ZERO_P (y))
  233.     ? (BIGNUM_MAYBE_COPY (x))
  234.     : ((BIGNUM_NEGATIVE_P (x))
  235.        ? ((BIGNUM_NEGATIVE_P (y))
  236.           ? (bignum_subtract_unsigned (y, x))
  237.           : (bignum_add_unsigned (x, y, 1)))
  238.        : ((BIGNUM_NEGATIVE_P (y))
  239.           ? (bignum_add_unsigned (x, y, 0))
  240.           : (bignum_subtract_unsigned (x, y))))));
  241. }
  242.  
  243. bignum_type
  244. DEFUN (bignum_negate, (x), fast bignum_type x)
  245. {
  246.   return
  247.     ((BIGNUM_ZERO_P (x))
  248.      ? (BIGNUM_MAYBE_COPY (x))
  249.      : (bignum_new_sign (x, (! (BIGNUM_NEGATIVE_P (x))))));
  250. }
  251.  
  252. bignum_type
  253. DEFUN (bignum_multiply, (x, y),
  254.        fast bignum_type x AND fast bignum_type y)
  255. {
  256.   fast bignum_length_type x_length = (BIGNUM_LENGTH (x));
  257.   fast bignum_length_type y_length = (BIGNUM_LENGTH (y));
  258.   fast int negative_p =
  259.     ((BIGNUM_NEGATIVE_P (x))
  260.      ? (! (BIGNUM_NEGATIVE_P (y)))
  261.      : (BIGNUM_NEGATIVE_P (y)));
  262.   if (BIGNUM_ZERO_P (x))
  263.     return (BIGNUM_MAYBE_COPY (x));
  264.   if (BIGNUM_ZERO_P (y))
  265.     return (BIGNUM_MAYBE_COPY (y));
  266.   if (x_length == 1)
  267.     {
  268.       bignum_digit_type digit = (BIGNUM_REF (x, 0));
  269.       if (digit == 1)
  270.     return (bignum_maybe_new_sign (y, negative_p));
  271.       if (digit < BIGNUM_RADIX_ROOT)
  272.     return (bignum_multiply_unsigned_small_factor (y, digit, negative_p));
  273.     }
  274.   if (y_length == 1)
  275.     {
  276.       bignum_digit_type digit = (BIGNUM_REF (y, 0));
  277.       if (digit == 1)
  278.     return (bignum_maybe_new_sign (x, negative_p));
  279.       if (digit < BIGNUM_RADIX_ROOT)
  280.     return (bignum_multiply_unsigned_small_factor (x, digit, negative_p));
  281.     }
  282.   return (bignum_multiply_unsigned (x, y, negative_p));
  283. }
  284.  
  285. int
  286. DEFUN (bignum_divide, (numerator, denominator, quotient, remainder),
  287.        bignum_type numerator AND bignum_type denominator
  288.        AND bignum_type * quotient AND bignum_type * remainder)
  289. {
  290.   if (BIGNUM_ZERO_P (denominator))
  291.     return (1);
  292.   if (BIGNUM_ZERO_P (numerator))
  293.     {
  294.       (*quotient) = (BIGNUM_MAYBE_COPY (numerator));
  295.       (*remainder) = (BIGNUM_MAYBE_COPY (numerator));
  296.     }
  297.   else
  298.     {
  299.       int r_negative_p = (BIGNUM_NEGATIVE_P (numerator));
  300.       int q_negative_p =
  301.     ((BIGNUM_NEGATIVE_P (denominator)) ? (! r_negative_p) : r_negative_p);
  302.       switch (bignum_compare_unsigned (numerator, denominator))
  303.     {
  304.     case bignum_comparison_equal:
  305.       {
  306.         (*quotient) = (BIGNUM_ONE (q_negative_p));
  307.         (*remainder) = (BIGNUM_ZERO ());
  308.         break;
  309.       }
  310.     case bignum_comparison_less:
  311.       {
  312.         (*quotient) = (BIGNUM_ZERO ());
  313.         (*remainder) = (BIGNUM_MAYBE_COPY (numerator));
  314.         break;
  315.       }
  316.     case bignum_comparison_greater:
  317.       {
  318.         if ((BIGNUM_LENGTH (denominator)) == 1)
  319.           {
  320.         bignum_digit_type digit = (BIGNUM_REF (denominator, 0));
  321.         if (digit == 1)
  322.           {
  323.             (*quotient) =
  324.               (bignum_maybe_new_sign (numerator, q_negative_p));
  325.             (*remainder) = (BIGNUM_ZERO ());
  326.             break;
  327.           }
  328.         else if (digit < BIGNUM_RADIX_ROOT)
  329.           {
  330.             bignum_divide_unsigned_small_denominator
  331.               (numerator, digit,
  332.                quotient, remainder,
  333.                q_negative_p, r_negative_p);
  334.             break;
  335.           }
  336.         else
  337.           {
  338.             bignum_divide_unsigned_medium_denominator
  339.               (numerator, digit,
  340.                quotient, remainder,
  341.                q_negative_p, r_negative_p);
  342.             break;
  343.           }
  344.           }
  345.         bignum_divide_unsigned_large_denominator
  346.           (numerator, denominator,
  347.            quotient, remainder,
  348.            q_negative_p, r_negative_p);
  349.         break;
  350.       }
  351.     }
  352.     }
  353.   return (0);
  354. }
  355.  
  356. bignum_type
  357. DEFUN (bignum_quotient, (numerator, denominator),
  358.        bignum_type numerator AND bignum_type denominator)
  359. {
  360.   if (BIGNUM_ZERO_P (denominator))
  361.     return (BIGNUM_OUT_OF_BAND);
  362.   if (BIGNUM_ZERO_P (numerator))
  363.     return (BIGNUM_MAYBE_COPY (numerator));
  364.   {
  365.     int q_negative_p =
  366.       ((BIGNUM_NEGATIVE_P (denominator))
  367.        ? (! (BIGNUM_NEGATIVE_P (numerator)))
  368.        : (BIGNUM_NEGATIVE_P (numerator)));
  369.     switch (bignum_compare_unsigned (numerator, denominator))
  370.       {
  371.       case bignum_comparison_equal:
  372.     return (BIGNUM_ONE (q_negative_p));
  373.       case bignum_comparison_less:
  374.     return (BIGNUM_ZERO ());
  375.       case bignum_comparison_greater:
  376.     {
  377.       bignum_type quotient;
  378.       if ((BIGNUM_LENGTH (denominator)) == 1)
  379.         {
  380.           bignum_digit_type digit = (BIGNUM_REF (denominator, 0));
  381.           if (digit == 1)
  382.         return (bignum_maybe_new_sign (numerator, q_negative_p));
  383.           if (digit < BIGNUM_RADIX_ROOT)
  384.         bignum_divide_unsigned_small_denominator
  385.           (numerator, digit,
  386.            ("ient), ((bignum_type *) 0),
  387.            q_negative_p, 0);
  388.           else
  389.         bignum_divide_unsigned_medium_denominator
  390.           (numerator, digit,
  391.            ("ient), ((bignum_type *) 0),
  392.            q_negative_p, 0);
  393.         }
  394.       else
  395.         bignum_divide_unsigned_large_denominator
  396.           (numerator, denominator,
  397.            ("ient), ((bignum_type *) 0),
  398.            q_negative_p, 0);
  399.       return (quotient);
  400.     }
  401.       }
  402.   }
  403. }
  404.  
  405. bignum_type
  406. DEFUN (bignum_remainder, (numerator, denominator),
  407.        bignum_type numerator AND bignum_type denominator)
  408. {
  409.   if (BIGNUM_ZERO_P (denominator))
  410.     return (BIGNUM_OUT_OF_BAND);
  411.   if (BIGNUM_ZERO_P (numerator))
  412.     return (BIGNUM_MAYBE_COPY (numerator));
  413.   switch (bignum_compare_unsigned (numerator, denominator))
  414.     {
  415.     case bignum_comparison_equal:
  416.       return (BIGNUM_ZERO ());
  417.     case bignum_comparison_less:
  418.       return (BIGNUM_MAYBE_COPY (numerator));
  419.     case bignum_comparison_greater:
  420.       {
  421.     bignum_type remainder;
  422.     if ((BIGNUM_LENGTH (denominator)) == 1)
  423.       {
  424.         bignum_digit_type digit = (BIGNUM_REF (denominator, 0));
  425.         if (digit == 1)
  426.           return (BIGNUM_ZERO ());
  427.         if (digit < BIGNUM_RADIX_ROOT)
  428.           return
  429.         (bignum_remainder_unsigned_small_denominator
  430.          (numerator, digit, (BIGNUM_NEGATIVE_P (numerator))));
  431.         bignum_divide_unsigned_medium_denominator
  432.           (numerator, digit,
  433.            ((bignum_type *) 0), (&remainder),
  434.            0, (BIGNUM_NEGATIVE_P (numerator)));
  435.       }
  436.     else
  437.       bignum_divide_unsigned_large_denominator
  438.         (numerator, denominator,
  439.          ((bignum_type *) 0), (&remainder),
  440.          0, (BIGNUM_NEGATIVE_P (numerator)));
  441.     return (remainder);
  442.       }
  443.     }
  444. }
  445.  
  446. /* These procedures depend on the non-portable type `unsigned long'.
  447.    If your compiler doesn't support this type, either define the
  448.    switch `BIGNUM_NO_ULONG' to disable them (in "bignum.h"), or write
  449.    new versions that don't use this type. */
  450.  
  451. #ifndef BIGNUM_NO_ULONG
  452.  
  453. bignum_type
  454. DEFUN (long_to_bignum, (n), long n)
  455. {
  456.   int negative_p;
  457.   bignum_digit_type result_digits [BIGNUM_DIGITS_FOR_LONG];
  458.   fast bignum_digit_type * end_digits = result_digits;
  459.   /* Special cases win when these small constants are cached. */
  460.   if (n == 0) return (BIGNUM_ZERO ());
  461.   if (n == 1) return (BIGNUM_ONE (0));
  462.   if (n == -1) return (BIGNUM_ONE (1));
  463.   {
  464.     fast unsigned long accumulator = ((negative_p = (n < 0)) ? (-n) : n);
  465.     do
  466.       {
  467.     (*end_digits++) = (accumulator & BIGNUM_DIGIT_MASK);
  468.     accumulator >>= BIGNUM_DIGIT_LENGTH;
  469.       }
  470.     while (accumulator != 0);
  471.   }
  472.   {
  473.     bignum_type result =
  474.       (bignum_allocate ((end_digits - result_digits), negative_p));
  475.     fast bignum_digit_type * scan_digits = result_digits;
  476.     fast bignum_digit_type * scan_result = (BIGNUM_START_PTR (result));
  477.     while (scan_digits < end_digits)
  478.       (*scan_result++) = (*scan_digits++);
  479.     return (result);
  480.   }
  481. }
  482.  
  483. long
  484. DEFUN (bignum_to_long, (bignum), bignum_type bignum)
  485. {
  486.   if (BIGNUM_ZERO_P (bignum))
  487.     return (0);
  488.   {
  489.     fast unsigned long accumulator = 0;
  490.     fast bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
  491.     fast bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum)));
  492.     while (start < scan)
  493.       accumulator = ((accumulator << BIGNUM_DIGIT_LENGTH) + (*--scan));
  494.     return ((BIGNUM_NEGATIVE_P (bignum)) ? (-accumulator) : accumulator);
  495.   }
  496. }
  497.  
  498. #endif /* not BIGNUM_NO_ULONG */
  499.  
  500. #define DTB_WRITE_DIGIT(factor)                        \
  501. {                                    \
  502.   significand *= (factor);                        \
  503.   digit = ((bignum_digit_type) significand);                \
  504.   (*--scan) = digit;                            \
  505.   significand -= ((double) digit);                    \
  506. }
  507.  
  508. bignum_type
  509. DEFUN (double_to_bignum, (x), double x)
  510. {
  511.   extern double frexp ();
  512.   int exponent;
  513.   fast double significand = (frexp (x, (&exponent)));
  514.   if (exponent <= 0) return (BIGNUM_ZERO ());
  515.   if (exponent == 1) return (BIGNUM_ONE (x < 0));
  516.   if (significand < 0) significand = (-significand);
  517.   {
  518.     bignum_length_type length = (BIGNUM_BITS_TO_DIGITS (exponent));
  519.     bignum_type result = (bignum_allocate (length, (x < 0)));
  520.     bignum_digit_type * start = (BIGNUM_START_PTR (result));
  521.     fast bignum_digit_type * scan = (start + length);
  522.     fast bignum_digit_type digit;
  523.     int odd_bits = (exponent % BIGNUM_DIGIT_LENGTH);
  524.     if (odd_bits > 0)
  525.       DTB_WRITE_DIGIT (1 << odd_bits);
  526.     while (start < scan)
  527.       {
  528.     if (significand == 0)
  529.       {
  530.         while (start < scan)
  531.           (*--scan) = 0;
  532.         break;
  533.       }
  534.     DTB_WRITE_DIGIT (BIGNUM_RADIX);
  535.       }
  536.     return (result);
  537.   }
  538. }
  539.  
  540. #undef DTB_WRITE_DIGIT
  541.  
  542. double
  543. DEFUN (bignum_to_double, (bignum), bignum_type bignum)
  544. {
  545.   if (BIGNUM_ZERO_P (bignum))
  546.     return (0);
  547.   {
  548.     fast double accumulator = 0;
  549.     fast bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
  550.     fast bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum)));
  551.     while (start < scan)
  552.       accumulator = ((accumulator * BIGNUM_RADIX) + (*--scan));
  553.     return ((BIGNUM_NEGATIVE_P (bignum)) ? (-accumulator) : accumulator);
  554.   }
  555. }
  556.  
  557. int
  558. DEFUN (bignum_fits_in_word_p, (bignum, word_length, twos_complement_p),
  559.        bignum_type bignum AND long word_length AND int twos_complement_p)
  560. {
  561.   unsigned int n_bits = (twos_complement_p ? (word_length - 1) : word_length);
  562.   BIGNUM_ASSERT (n_bits > 0);
  563.   {
  564.     fast bignum_length_type length = (BIGNUM_LENGTH (bignum));
  565.     fast unsigned int max_digits = (BIGNUM_BITS_TO_DIGITS (n_bits));
  566.     bignum_digit_type msd, max;
  567.     return
  568.       ((length < max_digits) ||
  569.        ((length == max_digits) &&
  570.     ((((msd = (BIGNUM_REF (bignum, (length - 1)))) <
  571.        (max = (1 << (n_bits - ((length - 1) * BIGNUM_DIGIT_LENGTH))))) ||
  572.       (twos_complement_p &&
  573.        (msd == max) &&
  574.        (BIGNUM_NEGATIVE_P (bignum)))))));
  575.   }
  576. }
  577.  
  578. bignum_type
  579. DEFUN (bignum_length_in_bits, (bignum), bignum_type bignum)
  580. {
  581.   if (BIGNUM_ZERO_P (bignum))
  582.     return (BIGNUM_ZERO ());
  583.   {
  584.     bignum_length_type index = ((BIGNUM_LENGTH (bignum)) - 1);
  585.     fast bignum_digit_type digit = (BIGNUM_REF (bignum, index));
  586.     fast bignum_type result = (bignum_allocate (2, 0));
  587.     (BIGNUM_REF (result, 0)) = index;
  588.     (BIGNUM_REF (result, 1)) = 0;
  589.     bignum_destructive_scale_up (result, BIGNUM_DIGIT_LENGTH);
  590.     while (digit > 0)
  591.       {
  592.     bignum_destructive_add (result, ((bignum_digit_type) 1));
  593.     digit >>= 1;
  594.       }
  595.     return (bignum_trim (result));
  596.   }
  597. }
  598.  
  599. bignum_type
  600. DEFUN_VOID (bignum_length_upper_limit)
  601. {
  602.   fast bignum_type result = (bignum_allocate (2, 0));
  603.   (BIGNUM_REF (result, 0)) = 0;
  604.   (BIGNUM_REF (result, 1)) = BIGNUM_DIGIT_LENGTH;
  605.   return (result);
  606. }
  607.  
  608. bignum_type
  609. DEFUN (digit_stream_to_bignum,
  610.        (n_digits, producer, context, radix, negative_p),
  611.        fast unsigned int n_digits
  612.        AND unsigned int EXFUN ((*producer), (bignum_procedure_context))
  613.        AND bignum_procedure_context context
  614.        AND fast unsigned int radix
  615.        AND int negative_p)
  616. {
  617.   BIGNUM_ASSERT ((radix > 1) && (radix <= BIGNUM_RADIX_ROOT));
  618.   if (n_digits == 0)
  619.     return (BIGNUM_ZERO ());
  620.   if (n_digits == 1)
  621.     {
  622.       long digit = ((long) ((*producer) (context)));
  623.       return (long_to_bignum (negative_p ? (- digit) : digit));
  624.     }
  625.   {
  626.     bignum_length_type length;
  627.     {
  628.       fast unsigned int radix_copy = radix;
  629.       fast unsigned int log_radix = 0;
  630.       while (radix_copy > 0)
  631.     {
  632.       radix_copy >>= 1;
  633.       log_radix += 1;
  634.     }
  635.       /* This length will be at least as large as needed. */
  636.       length = (BIGNUM_BITS_TO_DIGITS (n_digits * log_radix));
  637.     }
  638.     {
  639.       fast bignum_type result = (bignum_allocate_zeroed (length, negative_p));
  640.       while ((n_digits--) > 0)
  641.     {
  642.       bignum_destructive_scale_up (result, ((bignum_digit_type) radix));
  643.       bignum_destructive_add
  644.         (result, ((bignum_digit_type) ((*producer) (context))));
  645.     }
  646.       return (bignum_trim (result));
  647.     }
  648.   }
  649. }
  650.  
  651. void
  652. DEFUN (bignum_to_digit_stream, (bignum, radix, consumer, context),
  653.        bignum_type bignum
  654.        AND unsigned int radix
  655.        AND void EXFUN ((*consumer),
  656.                (bignum_procedure_context, bignum_digit_type))
  657.        AND bignum_procedure_context context)
  658. {
  659.   BIGNUM_ASSERT ((radix > 1) && (radix <= BIGNUM_RADIX_ROOT));
  660.   if (! (BIGNUM_ZERO_P (bignum)))
  661.     {
  662.       fast bignum_type working_copy = (bignum_copy (bignum));
  663.       fast bignum_digit_type * start = (BIGNUM_START_PTR (working_copy));
  664.       fast bignum_digit_type * scan = (start + (BIGNUM_LENGTH (working_copy)));
  665.       while (start < scan)
  666.     {
  667.       if ((scan[-1]) == 0)
  668.         scan -= 1;
  669.       else
  670.         (*consumer)
  671.           (context, (bignum_destructive_scale_down (working_copy, radix)));
  672.     }
  673.       BIGNUM_DEALLOCATE (working_copy);
  674.     }
  675.   return;
  676. }
  677.  
  678. long
  679. DEFUN_VOID (bignum_max_digit_stream_radix)
  680. {
  681.   return (BIGNUM_RADIX_ROOT);
  682. }
  683.  
  684. /* Comparisons */
  685.  
  686. static int
  687. DEFUN (bignum_equal_p_unsigned, (x, y),
  688.        bignum_type x AND bignum_type y)
  689. {
  690.   bignum_length_type length = (BIGNUM_LENGTH (x));
  691.   if (length != (BIGNUM_LENGTH (y)))
  692.     return (0);
  693.   else
  694.     {
  695.       fast bignum_digit_type * scan_x = (BIGNUM_START_PTR (x));
  696.       fast bignum_digit_type * scan_y = (BIGNUM_START_PTR (y));
  697.       fast bignum_digit_type * end_x = (scan_x + length);
  698.       while (scan_x < end_x)
  699.     if ((*scan_x++) != (*scan_y++))
  700.       return (0);
  701.       return (1);
  702.     }
  703. }
  704.  
  705. static enum bignum_comparison
  706. DEFUN (bignum_compare_unsigned, (x, y),
  707.        bignum_type x AND bignum_type y)
  708. {
  709.   bignum_length_type x_length = (BIGNUM_LENGTH (x));
  710.   bignum_length_type y_length = (BIGNUM_LENGTH (y));
  711.   if (x_length < y_length)
  712.     return (bignum_comparison_less);
  713.   if (x_length > y_length)
  714.     return (bignum_comparison_greater);
  715.   {
  716.     fast bignum_digit_type * start_x = (BIGNUM_START_PTR (x));
  717.     fast bignum_digit_type * scan_x = (start_x + x_length);
  718.     fast bignum_digit_type * scan_y = ((BIGNUM_START_PTR (y)) + y_length);
  719.     while (start_x < scan_x)
  720.       {
  721.     fast bignum_digit_type digit_x = (*--scan_x);
  722.     fast bignum_digit_type digit_y = (*--scan_y);
  723.     if (digit_x < digit_y)
  724.       return (bignum_comparison_less);
  725.     if (digit_x > digit_y)
  726.       return (bignum_comparison_greater);
  727.       }
  728.   }
  729.   return (bignum_comparison_equal);
  730. }
  731.  
  732. /* Addition */
  733.  
  734. static bignum_type
  735. DEFUN (bignum_add_unsigned, (x, y, negative_p),
  736.        bignum_type x AND bignum_type y AND int negative_p)
  737. {
  738.   if ((BIGNUM_LENGTH (y)) > (BIGNUM_LENGTH (x)))
  739.     {
  740.       bignum_type z = x;
  741.       x = y;
  742.       y = z;
  743.     }
  744.   {
  745.     bignum_length_type x_length = (BIGNUM_LENGTH (x));
  746.     bignum_type r = (bignum_allocate ((x_length + 1), negative_p));
  747.     fast bignum_digit_type sum;
  748.     fast bignum_digit_type carry = 0;
  749.     fast bignum_digit_type * scan_x = (BIGNUM_START_PTR (x));
  750.     fast bignum_digit_type * scan_r = (BIGNUM_START_PTR (r));
  751.     {
  752.       fast bignum_digit_type * scan_y = (BIGNUM_START_PTR (y));
  753.       fast bignum_digit_type * end_y = (scan_y + (BIGNUM_LENGTH (y)));
  754.       while (scan_y < end_y)
  755.     {
  756.       sum = ((*scan_x++) + (*scan_y++) + carry);
  757.       if (sum < BIGNUM_RADIX)
  758.         {
  759.           (*scan_r++) = sum;
  760.           carry = 0;
  761.         }
  762.       else
  763.         {
  764.           (*scan_r++) = (sum - BIGNUM_RADIX);
  765.           carry = 1;
  766.         }
  767.     }
  768.     }
  769.     {
  770.       fast bignum_digit_type * end_x = ((BIGNUM_START_PTR (x)) + x_length);
  771.       if (carry != 0)
  772.     while (scan_x < end_x)
  773.       {
  774.         sum = ((*scan_x++) + 1);
  775.         if (sum < BIGNUM_RADIX)
  776.           {
  777.         (*scan_r++) = sum;
  778.         carry = 0;
  779.         break;
  780.           }
  781.         else
  782.           (*scan_r++) = (sum - BIGNUM_RADIX);
  783.       }
  784.       while (scan_x < end_x)
  785.     (*scan_r++) = (*scan_x++);
  786.     }
  787.     if (carry != 0)
  788.       {
  789.     (*scan_r) = 1;
  790.     return (r);
  791.       }
  792.     return (bignum_shorten_length (r, x_length));
  793.   }
  794. }
  795.  
  796. /* Subtraction */
  797.  
  798. static bignum_type
  799. DEFUN (bignum_subtract_unsigned, (x, y),
  800.        bignum_type x AND bignum_type y)
  801. {
  802.   int negative_p;
  803.   switch (bignum_compare_unsigned (x, y))
  804.     {
  805.     case bignum_comparison_equal:
  806.       return (BIGNUM_ZERO ());
  807.     case bignum_comparison_less:
  808.       {
  809.     bignum_type z = x;
  810.     x = y;
  811.     y = z;
  812.       }
  813.       negative_p = 1;
  814.       break;
  815.     case bignum_comparison_greater:
  816.       negative_p = 0;
  817.       break;
  818.     }
  819.   {
  820.     bignum_length_type x_length = (BIGNUM_LENGTH (x));
  821.     bignum_type r = (bignum_allocate (x_length, negative_p));
  822.     fast bignum_digit_type difference;
  823.     fast bignum_digit_type borrow = 0;
  824.     fast bignum_digit_type * scan_x = (BIGNUM_START_PTR (x));
  825.     fast bignum_digit_type * scan_r = (BIGNUM_START_PTR (r));
  826.     {
  827.       fast bignum_digit_type * scan_y = (BIGNUM_START_PTR (y));
  828.       fast bignum_digit_type * end_y = (scan_y + (BIGNUM_LENGTH (y)));
  829.       while (scan_y < end_y)
  830.     {
  831.       difference = (((*scan_x++) - (*scan_y++)) - borrow);
  832.       if (difference < 0)
  833.         {
  834.           (*scan_r++) = (difference + BIGNUM_RADIX);
  835.           borrow = 1;
  836.         }
  837.       else
  838.         {
  839.           (*scan_r++) = difference;
  840.           borrow = 0;
  841.         }
  842.     }
  843.     }
  844.     {
  845.       fast bignum_digit_type * end_x = (scan_x + x_length);
  846.       if (borrow != 0)
  847.     while (scan_x < end_x)
  848.       {
  849.         difference = ((*scan_x++) - borrow);
  850.         if (difference < 0)
  851.           (*scan_r++) = (difference + BIGNUM_RADIX);
  852.         else
  853.           {
  854.         (*scan_r++) = difference;
  855.         borrow = 0;
  856.         break;
  857.           }
  858.       }
  859.       BIGNUM_ASSERT (borrow == 0);
  860.       while (scan_x < end_x)
  861.     (*scan_r++) = (*scan_x++);
  862.     }
  863.     return (bignum_trim (r));
  864.   }
  865. }
  866.  
  867. /* Multiplication
  868.    Maximum value for product_low or product_high:
  869.     ((R * R) + (R * (R - 2)) + (R - 1))
  870.    Maximum value for carry: ((R * (R - 1)) + (R - 1))
  871.     where R == BIGNUM_RADIX_ROOT */
  872.  
  873. static bignum_type
  874. DEFUN (bignum_multiply_unsigned, (x, y, negative_p),
  875.        bignum_type x AND bignum_type y AND int negative_p)
  876. {
  877.   if ((BIGNUM_LENGTH (y)) > (BIGNUM_LENGTH (x)))
  878.     {
  879.       bignum_type z = x;
  880.       x = y;
  881.       y = z;
  882.     }
  883.   {
  884.     fast bignum_digit_type carry;
  885.     fast bignum_digit_type y_digit_low;
  886.     fast bignum_digit_type y_digit_high;
  887.     fast bignum_digit_type x_digit_low;
  888.     fast bignum_digit_type x_digit_high;
  889.     bignum_digit_type product_low;
  890.     fast bignum_digit_type * scan_r;
  891.     fast bignum_digit_type * scan_y;
  892.     bignum_length_type x_length = (BIGNUM_LENGTH (x));
  893.     bignum_length_type y_length = (BIGNUM_LENGTH (y));
  894.     bignum_type r =
  895.       (bignum_allocate_zeroed ((x_length + y_length), negative_p));
  896.     bignum_digit_type * scan_x = (BIGNUM_START_PTR (x));
  897.     bignum_digit_type * end_x = (scan_x + x_length);
  898.     bignum_digit_type * start_y = (BIGNUM_START_PTR (y));
  899.     bignum_digit_type * end_y = (start_y + y_length);
  900.     bignum_digit_type * start_r = (BIGNUM_START_PTR (r));
  901. #define x_digit x_digit_high
  902. #define y_digit y_digit_high
  903. #define product_high carry
  904.     while (scan_x < end_x)
  905.       {
  906.     x_digit = (*scan_x++);
  907.     x_digit_low = (HD_LOW (x_digit));
  908.     x_digit_high = (HD_HIGH (x_digit));
  909.     carry = 0;
  910.     scan_y = start_y;
  911.     scan_r = (start_r++);
  912.     while (scan_y < end_y)
  913.       {
  914.         y_digit = (*scan_y++);
  915.         y_digit_low = (HD_LOW (y_digit));
  916.         y_digit_high = (HD_HIGH (y_digit));
  917.         product_low =
  918.           ((*scan_r) +
  919.            (x_digit_low * y_digit_low) +
  920.            (HD_LOW (carry)));
  921.         product_high =
  922.           ((x_digit_high * y_digit_low) +
  923.            (x_digit_low * y_digit_high) +
  924.            (HD_HIGH (product_low)) +
  925.            (HD_HIGH (carry)));
  926.         (*scan_r++) =
  927.           (HD_CONS ((HD_LOW (product_high)), (HD_LOW (product_low))));
  928.         carry =
  929.           ((x_digit_high * y_digit_high) +
  930.            (HD_HIGH (product_high)));
  931.       }
  932.     (*scan_r) += carry;
  933.       }
  934.     return (bignum_trim (r));
  935. #undef x_digit
  936. #undef y_digit
  937. #undef product_high
  938.   }
  939. }
  940.  
  941. static bignum_type
  942. DEFUN (bignum_multiply_unsigned_small_factor, (x, y, negative_p),
  943.        bignum_type x AND bignum_digit_type y AND int negative_p)
  944. {
  945.   bignum_length_type length_x = (BIGNUM_LENGTH (x));
  946.   bignum_type p = (bignum_allocate ((length_x + 1), negative_p));
  947.   bignum_destructive_copy (x, p);
  948.   (BIGNUM_REF (p, length_x)) = 0;
  949.   bignum_destructive_scale_up (p, y);
  950.   return (bignum_trim (p));
  951. }
  952.  
  953. static void
  954. DEFUN (bignum_destructive_scale_up, (bignum, factor),
  955.        bignum_type bignum AND bignum_digit_type factor)
  956. {
  957.   fast bignum_digit_type carry = 0;
  958.   fast bignum_digit_type * scan = (BIGNUM_START_PTR (bignum));
  959.   fast bignum_digit_type two_digits;
  960.   fast bignum_digit_type product_low;
  961. #define product_high carry
  962.   bignum_digit_type * end = (scan + (BIGNUM_LENGTH (bignum)));
  963.   BIGNUM_ASSERT ((factor > 1) && (factor < BIGNUM_RADIX_ROOT));
  964.   while (scan < end)
  965.     {
  966.       two_digits = (*scan);
  967.       product_low = ((factor * (HD_LOW (two_digits))) + (HD_LOW (carry)));
  968.       product_high =
  969.     ((factor * (HD_HIGH (two_digits))) +
  970.      (HD_HIGH (product_low)) +
  971.      (HD_HIGH (carry)));
  972.       (*scan++) = (HD_CONS ((HD_LOW (product_high)), (HD_LOW (product_low))));
  973.       carry = (HD_HIGH (product_high));
  974.     }
  975.   /* A carry here would be an overflow, i.e. it would not fit.
  976.      Hopefully the callers allocate enough space that this will
  977.      never happen.
  978.    */
  979.   BIGNUM_ASSERT (carry == 0);
  980.   return;
  981. #undef product_high
  982. }
  983.  
  984. static void
  985. DEFUN (bignum_destructive_add, (bignum, n),
  986.        bignum_type bignum AND bignum_digit_type n)
  987. {
  988.   fast bignum_digit_type * scan = (BIGNUM_START_PTR (bignum));
  989.   fast bignum_digit_type digit;
  990.   digit = ((*scan) + n);
  991.   if (digit < BIGNUM_RADIX)
  992.     {
  993.       (*scan) = digit;
  994.       return;
  995.     }
  996.   (*scan++) = (digit - BIGNUM_RADIX);
  997.   while (1)
  998.     {
  999.       digit = ((*scan) + 1);
  1000.       if (digit < BIGNUM_RADIX)
  1001.     {
  1002.       (*scan) = digit;
  1003.       return;
  1004.     }
  1005.       (*scan++) = (digit - BIGNUM_RADIX);
  1006.     }
  1007. }
  1008.  
  1009. /* Division */
  1010.  
  1011. /* For help understanding this algorithm, see:
  1012.    Knuth, Donald E., "The Art of Computer Programming",
  1013.    volume 2, "Seminumerical Algorithms"
  1014.    section 4.3.1, "Multiple-Precision Arithmetic". */
  1015.  
  1016. static void
  1017. DEFUN (bignum_divide_unsigned_large_denominator, (numerator, denominator,
  1018.                           quotient, remainder,
  1019.                           q_negative_p, r_negative_p),
  1020.        bignum_type numerator
  1021.        AND bignum_type denominator
  1022.        AND bignum_type * quotient
  1023.        AND bignum_type * remainder
  1024.        AND int q_negative_p
  1025.        AND int r_negative_p)
  1026. {
  1027.   bignum_length_type length_n = ((BIGNUM_LENGTH (numerator)) + 1);
  1028.   bignum_length_type length_d = (BIGNUM_LENGTH (denominator));
  1029.   bignum_type q =
  1030.     ((quotient != ((bignum_type *) 0))
  1031.      ? (bignum_allocate ((length_n - length_d), q_negative_p))
  1032.      : BIGNUM_OUT_OF_BAND);
  1033.   bignum_type u = (bignum_allocate (length_n, r_negative_p));
  1034.   int shift = 0;
  1035.   BIGNUM_ASSERT (length_d > 1);
  1036.   {
  1037.     fast bignum_digit_type v1 = (BIGNUM_REF ((denominator), (length_d - 1)));
  1038.     while (v1 < (BIGNUM_RADIX / 2))
  1039.       {
  1040.     v1 <<= 1;
  1041.     shift += 1;
  1042.       }
  1043.   }
  1044.   if (shift == 0)
  1045.     {
  1046.       bignum_destructive_copy (numerator, u);
  1047.       (BIGNUM_REF (u, (length_n - 1))) = 0;
  1048.       bignum_divide_unsigned_normalized (u, denominator, q);
  1049.     }
  1050.   else
  1051.     {
  1052.       bignum_type v = (bignum_allocate (length_d, 0));
  1053.       bignum_destructive_normalization (numerator, u, shift);
  1054.       bignum_destructive_normalization (denominator, v, shift);
  1055.       bignum_divide_unsigned_normalized (u, v, q);
  1056.       BIGNUM_DEALLOCATE (v);
  1057.       if (remainder != ((bignum_type *) 0))
  1058.     bignum_destructive_unnormalization (u, shift);
  1059.     }
  1060.   if (quotient != ((bignum_type *) 0))
  1061.     (*quotient) = (bignum_trim (q));
  1062.   if (remainder != ((bignum_type *) 0))
  1063.     (*remainder) = (bignum_trim (u));
  1064.   else
  1065.     BIGNUM_DEALLOCATE (u);
  1066.   return;
  1067. }
  1068.  
  1069. static void
  1070. DEFUN (bignum_divide_unsigned_normalized, (u, v, q),
  1071.        bignum_type u AND bignum_type v AND bignum_type q)
  1072. {
  1073.   bignum_length_type u_length = (BIGNUM_LENGTH (u));
  1074.   bignum_length_type v_length = (BIGNUM_LENGTH (v));
  1075.   bignum_digit_type * u_start = (BIGNUM_START_PTR (u));
  1076.   bignum_digit_type * u_scan = (u_start + u_length);
  1077.   bignum_digit_type * u_scan_limit = (u_start + v_length);
  1078.   bignum_digit_type * u_scan_start = (u_scan - v_length);
  1079.   bignum_digit_type * v_start = (BIGNUM_START_PTR (v));
  1080.   bignum_digit_type * v_end = (v_start + v_length);
  1081.   bignum_digit_type * q_scan;
  1082.   bignum_digit_type v1 = (v_end[-1]);
  1083.   bignum_digit_type v2 = (v_end[-2]);
  1084.   fast bignum_digit_type ph;    /* high half of double-digit product */
  1085.   fast bignum_digit_type pl;    /* low half of double-digit product */
  1086.   fast bignum_digit_type guess;
  1087.   fast bignum_digit_type gh;    /* high half-digit of guess */
  1088.   fast bignum_digit_type ch;    /* high half of double-digit comparand */
  1089.   fast bignum_digit_type v2l = (HD_LOW (v2));
  1090.   fast bignum_digit_type v2h = (HD_HIGH (v2));
  1091.   fast bignum_digit_type cl;    /* low half of double-digit comparand */
  1092. #define gl ph            /* low half-digit of guess */
  1093. #define uj pl
  1094. #define qj ph
  1095.   bignum_digit_type gm;        /* memory loc for reference parameter */
  1096.   if (q != BIGNUM_OUT_OF_BAND)
  1097.     q_scan = ((BIGNUM_START_PTR (q)) + (BIGNUM_LENGTH (q)));
  1098.   while (u_scan_limit < u_scan)
  1099.     {
  1100.       uj = (*--u_scan);
  1101.       if (uj != v1)
  1102.     {
  1103.       /* comparand =
  1104.          (((((uj * BIGNUM_RADIX) + uj1) % v1) * BIGNUM_RADIX) + uj2);
  1105.          guess = (((uj * BIGNUM_RADIX) + uj1) / v1); */
  1106.       cl = (u_scan[-2]);
  1107.       ch = (bignum_digit_divide (uj, (u_scan[-1]), v1, (&gm)));
  1108.       guess = gm;
  1109.     }
  1110.       else
  1111.     {
  1112.       cl = (u_scan[-2]);
  1113.       ch = ((u_scan[-1]) + v1);
  1114.       guess = (BIGNUM_RADIX - 1);
  1115.     }
  1116.       while (1)
  1117.     {
  1118.       /* product = (guess * v2); */
  1119.       gl = (HD_LOW (guess));
  1120.       gh = (HD_HIGH (guess));
  1121.       pl = (v2l * gl);
  1122.       ph = ((v2l * gh) + (v2h * gl) + (HD_HIGH (pl)));
  1123.       pl = (HD_CONS ((HD_LOW (ph)), (HD_LOW (pl))));
  1124.       ph = ((v2h * gh) + (HD_HIGH (ph)));
  1125.       /* if (comparand >= product) */
  1126.       if ((ch > ph) || ((ch == ph) && (cl >= pl)))
  1127.         break;
  1128.       guess -= 1;
  1129.       /* comparand += (v1 << BIGNUM_DIGIT_LENGTH) */
  1130.       ch += v1;
  1131.       /* if (comparand >= (BIGNUM_RADIX * BIGNUM_RADIX)) */
  1132.       if (ch >= BIGNUM_RADIX)
  1133.         break;
  1134.     }
  1135.       qj = (bignum_divide_subtract (v_start, v_end, guess, (--u_scan_start)));
  1136.       if (q != BIGNUM_OUT_OF_BAND)
  1137.     (*--q_scan) = qj;
  1138.     }
  1139.   return;
  1140. #undef gl
  1141. #undef uj
  1142. #undef qj
  1143. }
  1144.  
  1145. static bignum_digit_type
  1146. DEFUN (bignum_divide_subtract, (v_start, v_end, guess, u_start),
  1147.        bignum_digit_type * v_start
  1148.        AND bignum_digit_type * v_end
  1149.        AND bignum_digit_type guess
  1150.        AND bignum_digit_type * u_start)
  1151. {
  1152.   bignum_digit_type * v_scan = v_start;
  1153.   bignum_digit_type * u_scan = u_start;
  1154.   fast bignum_digit_type carry = 0;
  1155.   if (guess == 0) return (0);
  1156.   {
  1157.     bignum_digit_type gl = (HD_LOW (guess));
  1158.     bignum_digit_type gh = (HD_HIGH (guess));
  1159.     fast bignum_digit_type v;
  1160.     fast bignum_digit_type pl;
  1161.     fast bignum_digit_type vl;
  1162. #define vh v
  1163. #define ph carry
  1164. #define diff pl
  1165.     while (v_scan < v_end)
  1166.       {
  1167.     v = (*v_scan++);
  1168.     vl = (HD_LOW (v));
  1169.     vh = (HD_HIGH (v));
  1170.     pl = ((vl * gl) + (HD_LOW (carry)));
  1171.     ph = ((vl * gh) + (vh * gl) + (HD_HIGH (pl)) + (HD_HIGH (carry)));
  1172.     diff = ((*u_scan) - (HD_CONS ((HD_LOW (ph)), (HD_LOW (pl)))));
  1173.     if (diff < 0)
  1174.       {
  1175.         (*u_scan++) = (diff + BIGNUM_RADIX);
  1176.         carry = ((vh * gh) + (HD_HIGH (ph)) + 1);
  1177.       }
  1178.     else
  1179.       {
  1180.         (*u_scan++) = diff;
  1181.         carry = ((vh * gh) + (HD_HIGH (ph)));
  1182.       }
  1183.       }
  1184.     if (carry == 0)
  1185.       return (guess);
  1186.     diff = ((*u_scan) - carry);
  1187.     if (diff < 0)
  1188.       (*u_scan) = (diff + BIGNUM_RADIX);
  1189.     else
  1190.       {
  1191.     (*u_scan) = diff;
  1192.     return (guess);
  1193.       }
  1194. #undef vh
  1195. #undef ph
  1196. #undef diff
  1197.   }
  1198.   /* Subtraction generated carry, implying guess is one too large.
  1199.      Add v back in to bring it back down. */
  1200.   v_scan = v_start;
  1201.   u_scan = u_start;
  1202.   carry = 0;
  1203.   while (v_scan < v_end)
  1204.     {
  1205.       bignum_digit_type sum = ((*v_scan++) + (*u_scan) + carry);
  1206.       if (sum < BIGNUM_RADIX)
  1207.     {
  1208.       (*u_scan++) = sum;
  1209.       carry = 0;
  1210.     }
  1211.       else
  1212.     {
  1213.       (*u_scan++) = (sum - BIGNUM_RADIX);
  1214.       carry = 1;
  1215.     }
  1216.     }
  1217.   if (carry == 1)
  1218.     {
  1219.       bignum_digit_type sum = ((*u_scan) + carry);
  1220.       (*u_scan) = ((sum < BIGNUM_RADIX) ? sum : (sum - BIGNUM_RADIX));
  1221.     }
  1222.   return (guess - 1);
  1223. }
  1224.  
  1225. static void
  1226. DEFUN (bignum_divide_unsigned_medium_denominator, (numerator, denominator,
  1227.                            quotient, remainder,
  1228.                            q_negative_p, r_negative_p),
  1229.        bignum_type numerator
  1230.        AND bignum_digit_type denominator
  1231.        AND bignum_type * quotient
  1232.        AND bignum_type * remainder
  1233.        AND int q_negative_p
  1234.        AND int r_negative_p)
  1235. {
  1236.   bignum_length_type length_n = (BIGNUM_LENGTH (numerator));
  1237.   bignum_length_type length_q;
  1238.   bignum_type q;
  1239.   int shift = 0;
  1240.   /* Because `bignum_digit_divide' requires a normalized denominator. */
  1241.   while (denominator < (BIGNUM_RADIX / 2))
  1242.     {
  1243.       denominator <<= 1;
  1244.       shift += 1;
  1245.     }
  1246.   if (shift == 0)
  1247.     {
  1248.       length_q = length_n;
  1249.       q = (bignum_allocate (length_q, q_negative_p));
  1250.       bignum_destructive_copy (numerator, q);
  1251.     }
  1252.   else
  1253.     {
  1254.       length_q = (length_n + 1);
  1255.       q = (bignum_allocate (length_q, q_negative_p));
  1256.       bignum_destructive_normalization (numerator, q, shift);
  1257.     }
  1258.   {
  1259.     fast bignum_digit_type r = 0;
  1260.     fast bignum_digit_type * start = (BIGNUM_START_PTR (q)); 
  1261.     fast bignum_digit_type * scan = (start + length_q);
  1262.     bignum_digit_type qj;
  1263.     if (quotient != ((bignum_type *) 0))
  1264.       {
  1265.     while (start < scan)
  1266.       {
  1267.         r = (bignum_digit_divide (r, (*--scan), denominator, (&qj)));
  1268.         (*scan) = qj;
  1269.       }
  1270.     (*quotient) = (bignum_trim (q));
  1271.       }
  1272.     else
  1273.       {
  1274.     while (start < scan)
  1275.       r = (bignum_digit_divide (r, (*--scan), denominator, (&qj)));
  1276.     BIGNUM_DEALLOCATE (q);
  1277.       }
  1278.     if (remainder != ((bignum_type *) 0))
  1279.       {
  1280.     if (shift != 0)
  1281.       r >>= shift;
  1282.     (*remainder) = (bignum_digit_to_bignum (r, r_negative_p));
  1283.       }
  1284.   }
  1285.   return;
  1286. }
  1287.  
  1288. static void
  1289. DEFUN (bignum_destructive_normalization, (source, target, shift_left),
  1290.        bignum_type source AND bignum_type target AND int shift_left)
  1291. {
  1292.   fast bignum_digit_type digit;
  1293.   fast bignum_digit_type * scan_source = (BIGNUM_START_PTR (source));
  1294.   fast bignum_digit_type carry = 0;
  1295.   fast bignum_digit_type * scan_target = (BIGNUM_START_PTR (target));
  1296.   bignum_digit_type * end_source = (scan_source + (BIGNUM_LENGTH (source)));
  1297.   bignum_digit_type * end_target = (scan_target + (BIGNUM_LENGTH (target)));
  1298.   int shift_right = (BIGNUM_DIGIT_LENGTH - shift_left);
  1299.   bignum_digit_type mask = ((1 << shift_right) - 1);
  1300.   while (scan_source < end_source)
  1301.     {
  1302.       digit = (*scan_source++);
  1303.       (*scan_target++) = (((digit & mask) << shift_left) | carry);
  1304.       carry = (digit >> shift_right);
  1305.     }
  1306.   if (scan_target < end_target)
  1307.     (*scan_target) = carry;
  1308.   else
  1309.     BIGNUM_ASSERT (carry == 0);
  1310.   return;
  1311. }
  1312.  
  1313. static void
  1314. DEFUN (bignum_destructive_unnormalization, (bignum, shift_right),
  1315.        bignum_type bignum AND int shift_right)
  1316. {
  1317.   bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
  1318.   fast bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum)));
  1319.   fast bignum_digit_type digit;
  1320.   fast bignum_digit_type carry = 0;
  1321.   int shift_left = (BIGNUM_DIGIT_LENGTH - shift_right);
  1322.   bignum_digit_type mask = ((1 << shift_right) - 1);
  1323.   while (start < scan)
  1324.     {
  1325.       digit = (*--scan);
  1326.       (*scan) = ((digit >> shift_right) | carry);
  1327.       carry = ((digit & mask) << shift_left);
  1328.     }
  1329.   BIGNUM_ASSERT (carry == 0);
  1330.   return;
  1331. }
  1332.  
  1333. /* This is a reduced version of the division algorithm, applied to the
  1334.    case of dividing two bignum digits by one bignum digit.  It is
  1335.    assumed that the numerator and denominator are normalized. */
  1336.  
  1337. #define BDD_STEP(qn, j)                            \
  1338. {                                    \
  1339.   uj = (u[j]);                                \
  1340.   if (uj != v1)                                \
  1341.     {                                    \
  1342.       uj_uj1 = (HD_CONS (uj, (u[j + 1])));                \
  1343.       guess = (uj_uj1 / v1);                        \
  1344.       comparand = (HD_CONS ((uj_uj1 % v1), (u[j + 2])));        \
  1345.     }                                    \
  1346.   else                                    \
  1347.     {                                    \
  1348.       guess = (BIGNUM_RADIX_ROOT - 1);                    \
  1349.       comparand = (HD_CONS (((u[j + 1]) + v1), (u[j + 2])));        \
  1350.     }                                    \
  1351.   while ((guess * v2) > comparand)                    \
  1352.     {                                    \
  1353.       guess -= 1;                            \
  1354.       comparand += (v1 << BIGNUM_HALF_DIGIT_LENGTH);            \
  1355.       if (comparand >= BIGNUM_RADIX)                    \
  1356.     break;                                \
  1357.     }                                    \
  1358.   qn = (bignum_digit_divide_subtract (v1, v2, guess, (&u[j])));        \
  1359. }
  1360.  
  1361. static bignum_digit_type
  1362. DEFUN (bignum_digit_divide, (uh, ul, v, q),
  1363.        bignum_digit_type uh AND bignum_digit_type ul
  1364.        AND bignum_digit_type v AND bignum_digit_type * q) /* return value */
  1365. {
  1366.   fast bignum_digit_type guess;
  1367.   fast bignum_digit_type comparand;
  1368.   fast bignum_digit_type v1 = (HD_HIGH (v));
  1369.   fast bignum_digit_type v2 = (HD_LOW (v));
  1370.   fast bignum_digit_type uj;
  1371.   fast bignum_digit_type uj_uj1;
  1372.   bignum_digit_type q1;
  1373.   bignum_digit_type q2;
  1374.   bignum_digit_type u [4];
  1375.   if (uh == 0)
  1376.     {
  1377.       if (ul < v)
  1378.     {
  1379.       (*q) = 0;
  1380.       return (ul);
  1381.     }
  1382.       else if (ul == v)
  1383.     {
  1384.       (*q) = 1;
  1385.       return (0);
  1386.     }
  1387.     }
  1388.   (u[0]) = (HD_HIGH (uh));
  1389.   (u[1]) = (HD_LOW (uh));
  1390.   (u[2]) = (HD_HIGH (ul));
  1391.   (u[3]) = (HD_LOW (ul));
  1392.   v1 = (HD_HIGH (v));
  1393.   v2 = (HD_LOW (v));
  1394.   BDD_STEP (q1, 0);
  1395.   BDD_STEP (q2, 1);
  1396.   (*q) = (HD_CONS (q1, q2));
  1397.   return (HD_CONS ((u[2]), (u[3])));
  1398. }
  1399.  
  1400. #undef BDD_STEP
  1401.  
  1402. #define BDDS_MULSUB(vn, un, carry_in)                    \
  1403. {                                    \
  1404.   product = ((vn * guess) + carry_in);                    \
  1405.   diff = (un - (HD_LOW (product)));                    \
  1406.   if (diff < 0)                                \
  1407.     {                                    \
  1408.       un = (diff + BIGNUM_RADIX_ROOT);                    \
  1409.       carry = ((HD_HIGH (product)) + 1);                \
  1410.     }                                    \
  1411.   else                                    \
  1412.     {                                    \
  1413.       un = diff;                            \
  1414.       carry = (HD_HIGH (product));                    \
  1415.     }                                    \
  1416. }
  1417.  
  1418. #define BDDS_ADD(vn, un, carry_in)                    \
  1419. {                                    \
  1420.   sum = (vn + un + carry_in);                        \
  1421.   if (sum < BIGNUM_RADIX_ROOT)                        \
  1422.     {                                    \
  1423.       un = sum;                                \
  1424.       carry = 0;                            \
  1425.     }                                    \
  1426.   else                                    \
  1427.     {                                    \
  1428.       un = (sum - BIGNUM_RADIX_ROOT);                    \
  1429.       carry = 1;                            \
  1430.     }                                    \
  1431. }
  1432.  
  1433. static bignum_digit_type
  1434. DEFUN (bignum_digit_divide_subtract, (v1, v2, guess, u),
  1435.        bignum_digit_type v1 AND bignum_digit_type v2
  1436.        AND bignum_digit_type guess AND bignum_digit_type * u)
  1437. {
  1438.   {
  1439.     fast bignum_digit_type product;
  1440.     fast bignum_digit_type diff;
  1441.     fast bignum_digit_type carry;
  1442.     BDDS_MULSUB (v2, (u[2]), 0);
  1443.     BDDS_MULSUB (v1, (u[1]), carry);
  1444.     if (carry == 0)
  1445.       return (guess);
  1446.     diff = ((u[0]) - carry);
  1447.     if (diff < 0)
  1448.       (u[0]) = (diff + BIGNUM_RADIX);
  1449.     else
  1450.       {
  1451.     (u[0]) = diff;
  1452.     return (guess);
  1453.       }
  1454.   }
  1455.   {
  1456.     fast bignum_digit_type sum;
  1457.     fast bignum_digit_type carry;
  1458.     BDDS_ADD(v2, (u[2]), 0);
  1459.     BDDS_ADD(v1, (u[1]), carry);
  1460.     if (carry == 1)
  1461.       (u[0]) += 1;
  1462.   }
  1463.   return (guess - 1);
  1464. }
  1465.  
  1466. #undef BDDS_MULSUB
  1467. #undef BDDS_ADD
  1468.  
  1469. static void
  1470. DEFUN (bignum_divide_unsigned_small_denominator, (numerator, denominator,
  1471.                           quotient, remainder,
  1472.                           q_negative_p, r_negative_p),
  1473.        bignum_type numerator
  1474.        AND bignum_digit_type denominator
  1475.        AND bignum_type * quotient
  1476.        AND bignum_type * remainder
  1477.        AND int q_negative_p
  1478.        AND int r_negative_p)
  1479. {
  1480.   bignum_type q = (bignum_new_sign (numerator, q_negative_p));
  1481.   bignum_digit_type r = (bignum_destructive_scale_down (q, denominator));
  1482.   (*quotient) = (bignum_trim (q));
  1483.   if (remainder != ((bignum_type *) 0))
  1484.     (*remainder) = (bignum_digit_to_bignum (r, r_negative_p));
  1485.   return;
  1486. }
  1487.  
  1488. /* Given (denominator > 1), it is fairly easy to show that
  1489.    (quotient_high < BIGNUM_RADIX_ROOT), after which it is easy to see
  1490.    that all digits are < BIGNUM_RADIX. */
  1491.  
  1492. static bignum_digit_type
  1493. DEFUN (bignum_destructive_scale_down, (bignum, denominator),
  1494.        bignum_type bignum AND fast bignum_digit_type denominator)
  1495. {
  1496.   fast bignum_digit_type numerator;
  1497.   fast bignum_digit_type remainder = 0;
  1498.   fast bignum_digit_type two_digits;
  1499. #define quotient_high remainder
  1500.   bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
  1501.   bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum)));
  1502.   BIGNUM_ASSERT ((denominator > 1) && (denominator < BIGNUM_RADIX_ROOT));
  1503.   while (start < scan)
  1504.     {
  1505.       two_digits = (*--scan);
  1506.       numerator = (HD_CONS (remainder, (HD_HIGH (two_digits))));
  1507.       quotient_high = (numerator / denominator);
  1508.       numerator = (HD_CONS ((numerator % denominator), (HD_LOW (two_digits))));
  1509.       (*scan) = (HD_CONS (quotient_high, (numerator / denominator)));
  1510.       remainder = (numerator % denominator);
  1511.     }
  1512.   return (remainder);
  1513. #undef quotient_high
  1514. }
  1515.  
  1516. static bignum_type
  1517. DEFUN (bignum_remainder_unsigned_small_denominator, (n, d, negative_p),
  1518.        bignum_type n AND bignum_digit_type d AND int negative_p)
  1519. {
  1520.   fast bignum_digit_type two_digits;
  1521.   bignum_digit_type * start = (BIGNUM_START_PTR (n));
  1522.   fast bignum_digit_type * scan = (start + (BIGNUM_LENGTH (n)));
  1523.   fast bignum_digit_type r = 0;
  1524.   BIGNUM_ASSERT ((d > 1) && (d < BIGNUM_RADIX_ROOT));
  1525.   while (start < scan)
  1526.     {
  1527.       two_digits = (*--scan);
  1528.       r =
  1529.     ((HD_CONS (((HD_CONS (r, (HD_HIGH (two_digits)))) % d),
  1530.            (HD_LOW (two_digits))))
  1531.      % d);
  1532.     }
  1533.   return (bignum_digit_to_bignum (r, negative_p));
  1534. }
  1535.  
  1536. static bignum_type
  1537. DEFUN (bignum_digit_to_bignum, (digit, negative_p),
  1538.        fast bignum_digit_type digit AND int negative_p)
  1539. {
  1540.   if (digit == 0)
  1541.     return (BIGNUM_ZERO ());
  1542.   else
  1543.     {
  1544.       fast bignum_type result = (bignum_allocate (1, negative_p));
  1545.       (BIGNUM_REF (result, 0)) = digit;
  1546.       return (result);
  1547.     }
  1548. }
  1549.  
  1550. /* Allocation */
  1551.  
  1552. static bignum_type
  1553. DEFUN (bignum_allocate, (length, negative_p),
  1554.        fast bignum_length_type length AND int negative_p)
  1555. {
  1556.   BIGNUM_ASSERT ((length >= 0) || (length < BIGNUM_RADIX));
  1557.   {
  1558.     fast bignum_type result = (BIGNUM_ALLOCATE (length));
  1559.     BIGNUM_SET_HEADER (result, length, negative_p);
  1560.     return (result);
  1561.   }
  1562. }
  1563.  
  1564. static bignum_type
  1565. DEFUN (bignum_allocate_zeroed, (length, negative_p),
  1566.        fast bignum_length_type length AND int negative_p)
  1567. {
  1568.   BIGNUM_ASSERT ((length >= 0) || (length < BIGNUM_RADIX));
  1569.   {
  1570.     fast bignum_type result = (BIGNUM_ALLOCATE (length));
  1571.     fast bignum_digit_type * scan = (BIGNUM_START_PTR (result));
  1572.     fast bignum_digit_type * end = (scan + length);
  1573.     BIGNUM_SET_HEADER (result, length, negative_p);
  1574.     while (scan < end)
  1575.       (*scan++) = 0;
  1576.     return (result);
  1577.   }
  1578. }
  1579.  
  1580. static bignum_type
  1581. DEFUN (bignum_shorten_length, (bignum, length),
  1582.        fast bignum_type bignum AND fast bignum_length_type length)
  1583. {
  1584.   fast bignum_length_type current_length = (BIGNUM_LENGTH (bignum));
  1585.   BIGNUM_ASSERT ((length >= 0) || (length <= current_length));
  1586.   if (length < current_length)
  1587.     {
  1588.       BIGNUM_SET_HEADER
  1589.     (bignum, length, ((length != 0) && (BIGNUM_NEGATIVE_P (bignum))));
  1590.       BIGNUM_REDUCE_LENGTH (bignum, bignum, length);
  1591.     }
  1592.   return (bignum);
  1593. }
  1594.  
  1595. static bignum_type
  1596. DEFUN (bignum_trim, (bignum), bignum_type bignum)
  1597. {
  1598.   fast bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
  1599.   fast bignum_digit_type * end = (start + (BIGNUM_LENGTH (bignum)));
  1600.   fast bignum_digit_type * scan = end;
  1601.   while ((start <= scan) && ((*--scan) == 0))
  1602.     ;
  1603.   scan += 1;
  1604.   if (scan < end)
  1605.     {
  1606.       fast bignum_length_type length = (scan - start);
  1607.       BIGNUM_SET_HEADER
  1608.     (bignum, length, ((length != 0) && (BIGNUM_NEGATIVE_P (bignum))));
  1609.       BIGNUM_REDUCE_LENGTH (bignum, bignum, length);
  1610.     }
  1611.   return (bignum);
  1612. }
  1613.  
  1614. /* Copying */
  1615.  
  1616. static bignum_type
  1617. DEFUN (bignum_copy, (source), fast bignum_type source)
  1618. {
  1619.   fast bignum_type target =
  1620.     (bignum_allocate ((BIGNUM_LENGTH (source)), (BIGNUM_NEGATIVE_P (source))));
  1621.   bignum_destructive_copy (source, target);
  1622.   return (target);
  1623. }
  1624.  
  1625. static bignum_type
  1626. DEFUN (bignum_new_sign, (bignum, negative_p),
  1627.        fast bignum_type bignum AND int negative_p)
  1628. {
  1629.   fast bignum_type result =
  1630.     (bignum_allocate ((BIGNUM_LENGTH (bignum)), negative_p));
  1631.   bignum_destructive_copy (bignum, result);
  1632.   return (result);
  1633. }
  1634.  
  1635. static bignum_type
  1636. DEFUN (bignum_maybe_new_sign, (bignum, negative_p),
  1637.        fast bignum_type bignum AND int negative_p)
  1638. {
  1639. #ifndef BIGNUM_FORCE_NEW_RESULTS
  1640.   if ((BIGNUM_NEGATIVE_P (bignum)) ? negative_p : (! negative_p))
  1641.     return (bignum);
  1642.   else
  1643. #endif /* not BIGNUM_FORCE_NEW_RESULTS */
  1644.     {
  1645.       fast bignum_type result =
  1646.     (bignum_allocate ((BIGNUM_LENGTH (bignum)), negative_p));
  1647.       bignum_destructive_copy (bignum, result);
  1648.       return (result);
  1649.     }
  1650. }
  1651.  
  1652. static void
  1653. DEFUN (bignum_destructive_copy, (source, target),
  1654.        bignum_type source AND bignum_type target)
  1655. {
  1656.   fast bignum_digit_type * scan_source = (BIGNUM_START_PTR (source));
  1657.   fast bignum_digit_type * end_source =
  1658.     (scan_source + (BIGNUM_LENGTH (source)));
  1659.   fast bignum_digit_type * scan_target = (BIGNUM_START_PTR (target));
  1660.   while (scan_source < end_source)
  1661.     (*scan_target++) = (*scan_source++);
  1662.   return;
  1663. }
  1664.  
  1665. static void
  1666. DEFUN (bignum_destructive_zero, (bignum), fast bignum_type bignum)
  1667. {
  1668.   fast bignum_digit_type * scan = (BIGNUM_START_PTR (bignum));
  1669.   fast bignum_digit_type * end = (scan + (BIGNUM_LENGTH (bignum)));
  1670.   while (scan < end)
  1671.     (*scan++) = 0;
  1672.   return;
  1673. }
  1674.